{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inference in Graphical Models\n",
    "\n",
    "The marginal inference of a graphical model is seeked to answer a question about a random variable of the graph. More specifically,\n",
    "\n",
    "$$\n",
    "    p(y\\in D) = \\sum_{x_1}\\sum_{x_2}\\cdots\\sum_{x_N}p(y\\in D, x_1, x_2, \\ldots, x_N)\n",
    "$$\n",
    "\n",
    "* Inference is often an NP-hard problem\n",
    "* Tractability of the inference is dependent on the structure of the graph that describs that probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "## Variable Elimination on a Chain\n",
    "Consider a Markov Random Field with discrete random variables $\\{x_n\\}_{n=1}^N$ taking $K$ elements each. The join probability of this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3wAAAJcCAYAAACiz6agAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3RUdf7/8ddMei+EUJLQRBAymZBg21VXXVfXFdxF5bfFVZJQxAqIIqACuioiPdIJkCJVuvTeu/SiQgQpoQkmgUBCyszvD7/LihQhJLkzk+fjHI/JnXs/n1c8R/E9r8m9JrvdbhcAAAAAwOWYjQ4AAAAAACgfDHwAAAAA4KIY+AAAAADARTHwAQAAAICLYuADAAAAABfFwAcAAAAALoqBDwAAAABcFAMfAAAAALgoBj4AAAAAcFEMfAAAAADgohj4AAAAAMBFMfABAAAAgIti4AMAAAAAF8XABwAAAAAuioEPAAAAAFwUAx8AAAAAuCgGPgAAAABwUQx8AAAAAOCiGPgAAAAAwEUx8AEAAACAi2LgAwAAAAAXxcAHAAAAAC6KgQ8AAAAAXBQDHwAAAAC4KAY+AAAAAHBRDHwAAAAA4KIY+AAAAADARTHwAQAAAICLcjc6QEU7mn1Ru4/landWrvYeP6fc/EIVltjl6WZSkI+nomsGKiYiSDGRQYoK8TU6LgAAAACUmslut9uNDlHeCopKNGfXcY3feFg7j+Xe9HWxkUF68f46am6tIW8Pt3JMCAAAAABlz6UHvhKbXePWHdKwFZnKyS8q9TrBPh567dH6av1AXbmZTWWYEAAAAADKj8sOfJmn89Rl2k5tP5pTZmvGRQWrX8tY1Q/3L7M1AQAAAKC8uOTAN2t7lt6ZvkuFJbYyX9vTzay+z1nVIi6izNcGAAAAgLLkcgPfxM1H9O7M3eW6h0nSJ8/E6Pl7a5XrPgAAAABwO1xq4Ju1PUudvtxRIXuZJA36exOaPgAAAAAOy2Wew5d5Ok/vTN9VYfvZJb0zfZcyT+dV2J4AAAAAcCtcouErsdnVcuT6Ut+gpTjnlLJGtvnN86onDJJXjTuvOBZfK1hT2/+eu3cCAAAAcDgu0fCNW3eozO7GafYLvu5fJvPVz+LbdiRH49YdKpO9AQAAAKAsuRsd4HYVFJVo2IrMMlsv6o3xt3zNsBWZevH+2jycHQAAAIBDcfqGb86u47f1UPWykJNfpLm7ThiaAQAAAAB+zekHvvEbD1/+Ontlmg73aa6jyc+rJC/7muefmtJTh/s014nUjrKXFJdZji9+kQMAAAAAHIFTD3xHsy9q57Hcy98HP/SCPKvdIVv+OZ2ZP/iq889tnaOCQ9tkcvdS2NNvy+RWdp9o3XksR0ezL5bZegAAAABwu5z6d/h2/2LYkySTm7vCnn5bJ9I6qeDgVp3fNk8B8c0kSUVnjylnRZokKfiRRHmERV1zzRMZb6nozBHJViKzX4i8IxvJv8lf5B0V/Zt59mTlKirE9/Z+KAAAAAAoI07d8O3Oyr3qmEdYlIIfSZQkZS8fp6Kzx2S3lejMnAGyF1+Sd504BTRtft01C49/J5l+/sdSkntKF/au1KkJXfXT0tH6rSdYXCsPAAAAABjFqRu+vcfPXfN4QNPmyv9+iwoObdOZOQPkXbeJCk8ekNk7QFWadZLJdOUz80zunvKPbya/Rg/Js9odMnv6yG63q/DU98pdO1H5mZt1/uuv5OYbrKDf//26efZkXTsPAAAAABjBqRu+3PzCax43mUyq0qyTzD6BKjx5QOc2TJUkhf75VbkHVLnqfDf/EFV54hV5R1lk9vS5vIZX9foKb9lTvnc9+PN+G76UrSDvBnmMvVsoAAAAAPySUw98hSXX/4ilu3+ogh9udfl737selF+jh0q1z38/ImovKlD+DztvkMdWqvUBAAAAoDw49cDn6Wa67mt2W4ku7Fl++fvCUwdlKywo1T4ewdVl9g2SJBXnnrxBHqf+xwkAAADAxTj1hBLk43nd185tnKZLx/bJ5OUnt8CqKs4+ruzlY0q/2eUbtlx/yAzy8Sj9+gAAAABQxpx64IuuGXjN45dOZipn7SRJUujj7RXW7E1JJuXtWKj877fc8j5FOSdly//5hizuQeHXPc8Sce08AAAAAGAEpx74YiKCrjpmK7qks3MGSLZi+TZ8QP6WP8q7tlUB9/xNknR2/ucquXjl4xN+63ELOasyJEkmdy951469pTwAAAAAYBTnHvgirx6wclamqejsUbn5hyr0ydcuHw95OEEeYbVUciFbZxcOveKaUxO7K3fDlyr88QfZbSWSfh4CC09+r9PTP9bFb1ZLkgLvf05uPgHXzWNh4AMAAADgQJz6OXxRIb6KjQzSzmM/N3b5h7br/Na5kqQqT3WUm8//PmJpcvdQ2NNv60RGZ+Xv36C8XUvlb/2TJKk497RyVmX83OSZ3WX28pW96JLsxZcuXx8Q31xBD/zrulliI4MVFeJbHj8mAAAAAJSKUzd8kvTC/bUlSSUFeTo7f7Aku/zjm8mnXtOrzvWsVk/BDz4vSfpp2WgV55ySJIX8sbX8Y/8sj/C6Mvv4y3bpomQyyT00Un7Wx1W91QCFPvHyVQ9s/6UX/y8HAAAAADgKk/23foHNwRUUlej+T5cpx8CHngf7eGhj98fk7eFmWAYAAAAA+DWnb/i8Pdz02qP1jdn8/2bl1x6tz7AHAAAAwOE4/cAnSa0fqKu4qOCK39hkkv3Hg6pbeKji9wYAAACA3+ASA5+b2aR+LWPl6VbBP05JsY7P/EyP/+lP6tChgy5evFix+wMAAADADbjEwCdJ9cP91fc5q65/W5WyZZLU/+9x6tHxJbm7u2vIkCFq0qSJNmzYUEEJAAAAAODGXGbgk6QWcRH65JmYch/6TJI+eSZGLe+upZ49e2rTpk2Kjo7WgQMH9OCDD6p79+66dOnSb64DAAAAAOXJ6e/SeS2ztmfpnem7VFhiK/O1vdzN+uxZq1rERVxx/NKlS+rVq5f69esnm80mi8WijIwMxcXFlXkGAAAAALgZLjnwSVLm6Tx1mbZT24/mlNma8bWC1a9lrO6o6n/dc9avX6+EhARlZmbK3d1dPXv2VLdu3eTh4VFmOQAAAADgZrjswCdJJTa7xq07pGErMm/rOX3BPh567dH6av1AXbmZf/sDoxcuXFD37t01ZMgQSdLdd9+t9PR0NW7cuNQZAAAAAOBWufTA918FRSWau+uEvth4WDuP3XzjFxsZrFa/q61mMTVK9Zy95cuXKykpSUeOHJGXl5d69+6tjh07ys2NZ/YBAAAAKH+VYuD7paPZF7UnK1e7s3K1J+uccvOLVFhik6ebWUE+HrJEBComIkiWiCBFhfje9n7nzp1T586dNXbsWEnSQw89pNTUVN1xxx23vTYAAAAA3EilG/iMMnfuXLVr104nT56Un5+f+vfvr/bt28tkqqgHSQAAAACobBj4KtDZs2f1+uuva/LkyZKkJ554QmPGjFFUVJTByQAAAAC4Ipd6Dp+jq1KliiZNmqQpU6aoSpUqWrx4sWJiYpSRkSHmbgAAAABljYbPICdPntRLL72kOXPmSJJatGihkSNHqlq1agYnAwAAAOAqaPgMUr16dc2ePVtpaWkKDAzUrFmzZLFYNG3aNKOjAQAAAHARNHwO4OjRo2rdurWWLl0qSXr++ec1ZMgQhYaGGpwMAAAAgDOj4XMAUVFRWrRokYYNGyZfX19NnDhRFotF8+fPNzoaAAAAACdGw+dgMjMzlZiYqHXr1kmS2rZtqwEDBigwMNDgZAAAAACcDQ2fg6lfv75WrVqlfv36ydPTU2PGjJHVatWKFSuMjgYAAADAydDwObC9e/eqVatW2rZtmySpQ4cO+vTTT+Xr62twMgAAAADOgIHPwRUVFal37976+OOPVVxcrAYNGig9PV3333+/0dEAAAAAODgGPiexdetWJSQkaO/evTKbzeratat69eolLy8vo6MBAAAAcFAMfE6koKBAvXr1Ur9+/WS32xUTE6OMjAw1adLE6GgAAAAAHBADnxNav369EhISlJmZKXd3d/Xq1UvdunWTu7u70dEAAAAAOBAGPid14cIFdevWTUOHDpUk3XPPPUpPT1ejRo0MTgYAAADAUfBYBifl5+enIUOGaOnSpYqKitKWLVsUFxengQMHymazGR0PAAAAgAOg4XMBubm5evPNN5WamipJ+sMf/qDU1FTVq1fP4GQAAAAAjETD5wKCgoI0btw4ffXVV6pWrZpWr14tq9WqUaNGiXkeAAAAqLxo+FzM2bNn9dprr2nKlCmSpCeeeEJjx45VZGSkwckAAAAAVDQaPhdTpUoVTZ48WZMnT1ZoaKgWL14si8WiL774grYPAAAAqGRo+FzYyZMn1a5dO82dO1eS1KJFC40aNUrh4eEGJwMAAABQEWj4XFj16tX11VdfKTU1VYGBgZo1a5aio6M1Y8YMo6MBAAAAqAA0fJXEkSNHlJSUpOXLl0uS/v3vf2vIkCEKCQkxOBkAAACA8kLDV0nUqlVLS5Ys0dChQ+Xr66sJEybIYrFowYIFRkcDAAAAUE5o+CqhAwcOKDExUevXr5cktWvXTgMGDFBAQIDByQAAAACUJRq+SujOO+/U6tWr1bdvX3l6eiolJUVWq1UrV640OhoAAACAMkTDV8nt2bNHCQkJ2rZtmySpY8eO+vTTT+Xj42NwMgAAAAC3i4avkrNYLNq4caN69eolNzc3JScnKy4uTps2bTI6GgAAAIDbRMOHy7Zu3apWrVpp3759MpvN6tatm3r27CkvLy+jowEAAAAoBQY+XKGgoEA9e/ZU//79ZbfbZbValZGRodjYWKOjAQAAALhFDHy4pnXr1ikhIUHff/+9PDw81KtXL3Xt2lXu7u5GRwMAAABwkxj4cF0XLlxQ165dNWzYMEnSPffco4yMDN11110GJwMAAABwM7hpC67Lz89PQ4cO1ZIlSxQVFaUtW7YoLi5OgwYNks1mMzoeAAAAgN9Aw4ebkpubq06dOiktLU2S9Ic//EFpaWmqW7euscEAAAAAXBcNH25KUFCQUlNTNXv2bFWrVk2rV69WTEyMRo8eLd4zAAAAABwTDR9u2ZkzZ/Tqq69q6tSpkqQ///nPGjt2rCIiIgxOBgAAAOCXaPhwy8LCwvTll19q8uTJCg0N1aJFi2SxWDR+/HjaPgAAAMCB0PDhtpw4cULt2rXTvHnzJEnPPvusRowYofDwcIOTAQAAAKDhw22pUaOG5syZo7FjxyogIEAzZsyQxWLRjBkzjI4GAAAAVHo0fCgzhw8fVuvWrbV8+XJJ0gsvvKDPP/9cISEhBicDAAAAKicaPpSZ2rVra8mSJRoyZIh8fHw0fvx4WSwWLVy40OhoAAAAQKVEw4dyceDAASUkJGjDhg2SpJdeekn9+/dXQECAwckAAACAyoOGD+Xizjvv1Jo1a/TZZ5/J09NTo0ePltVq1apVq4yOBgAAAFQaNHwod3v27FGrVq20fft2mUwmdezYUb1795aPj4/R0QAAAACXRsOHcmexWLRx40b17NlTZrNZgwcPVlxcnDZv3mx0NAAAAMCl0fChQn399ddq1aqVvvnmG5nNZnXv3l09e/aUp6en0dEAAAAAl8PAhwpXUFCgHj16aMCAAbLb7bJarcrIyFBsbKzR0QAAAACXwsAHw6xdu1YJCQk6ePCgPDw89MEHH+idd96Ru7u70dEAAAAAl8DAB0Pl5eWpa9euGj58uCTp3nvvVXp6uu666y6DkwEAAADOj5u2wFD+/v4aNmyYFi1apMjISG3evFlxcXEaPHiwbDab0fEAAAAAp0bDB4eRk5OjTp06KT09XZL08MMPKzU1VXXr1jU4GQAAAOCcaPjgMIKDg5WWlqZZs2YpPDxcq1atktVqVUpKinhfAgAAALh1NHxwSGfOnNErr7yiadOmSZKefPJJjRkzRhEREQYnAwAAAJwHDR8cUlhYmL788ktNnDhRISEhWrhwoSwWiyZMmEDbBwAAANwkGj44vOPHj6tdu3aaP3++JOnZZ5/VyJEjVbVqVYOTAQAAAI6Nhg8Or2bNmpo7d67GjBmjgIAAzZgxQ9HR0Zo1a5bR0QAAAACHRsMHp/LDDz+odevWWrFihSTpxRdfVHJyskJCQgxOBgAAADgeGj44lTp16mjp0qX6/PPP5ePjoy+++EIxMTFatGiR0dEAAAAAh0PDB6e1f/9+JSQkaOPGjZKk9u3bq1+/fgoICDA4GQAAAOAYaPjgtBo0aKC1a9eqT58+8vT01KhRoxQbG6vVq1cbHQ0AAABwCDR8cAm7d+9Wq1attGPHDplMJnXq1EmffPKJfHx8jI4GAAAAGIaGDy4hJiZGmzZtUo8ePWQ2mzVo0CDFx8dr8+bNRkcDAAAADEPDB5ezZcsWtWrVSt9++63c3NzUvXt39ejRQ56enkZHAwAAACoUAx9cUn5+vnr06KGBAwfKbrcrNjZWGRkZslqtRkcDAAAAKgwDH1zamjVrlJiYqIMHD8rDw0MffvihunTpInd3d6OjAQAAAOWOgQ8uLy8vT126dNHIkSMlSffdd5/S09PVsGFDg5MBAAAA5YubtsDl+fv7a8SIEVq4cKEiIiK0adMmxcXFKTk5WTabzeh4AAAAQLmh4UOlkpOTo44dOyojI0OS9Mgjjyg1NVV16tQxNhgAAABQDmj4UKkEBwcrPT1dM2fOVHh4uFauXKmYmBiNGTNGvPcBAAAAV0PDh0rrxx9/1CuvvKLp06dLkv7yl79ozJgxqlmzpsHJAAAAgLJBw4dKq2rVqpo6daomTJig4OBgLViwQBaLRRMnTqTtAwAAgEug4QMkHT9+XG3bttWCBQskSc8995xGjBihqlWrGpwMAAAAKD0aPkBSzZo1NW/ePKWkpMjf31/Tp09XdHS0Zs2aZXQ0AAAAoNRo+IBf+eGHH5SUlKSVK1dKklq1aqXk5GQFBwcbGwwAAAC4RTR8wK/UqVNHy5YtU3Jysry9vZWRkSGLxaLFixcbHQ0AAAC4JTR8wA189913SkhI0KZNmyRJL7/8svr16yd/f3+DkwEAAAC/jYYPuIGGDRtq7dq16t27tzw8PDRy5EjFxsZqzZo1RkcDAAAAfhMNH3CTdu3apVatWmnnzp0ymUzq3LmzPvroI/n4+BgdDQAAALgmGj7gJlmtVm3evFnvv/++zGazBgwYoKZNm2rLli1GRwMAAACuiYYPKIXNmzcrISFB3377rdzc3PTuu+/q/fffl6enp9HRAAAAgMsY+IBSys/P13vvvafBgwfLbrerSZMmysjIUExMjNHRAAAAAEl8pBMoNR8fHw0cOFArVqxQ3bp1tWPHDjVt2lR9+vRRSUmJ0fEAAAAAGj6gLJw/f15dunTRqFGjJEn333+/0tPT1aBBA4OTAQAAoDJj4APK0KJFi9SmTRtlZWXJx8dHffr00euvvy6z+fbL9KPZF7X7WK52Z+Vq7/Fzys0vVGGJXZ5uJgX5eCq6ZqBiIoIUExmkqBDfMvhpAAAA4OwY+IAylp2drY4dO+qLL76QJD366KNKTU1V7dq1b3mtgqISzdl1XOM3HtbOY7k3fV1sZJBevL+OmltryNvD7Zb3BQAAgGtg4APKycyZM9W+fXv9+OOPCggI0KBBg9S6dWuZTKbfvLbEZte4dYc0bEWmcvKLSp0h2MdDrz1aX60fqCs382/vCwAAANfCwAeUox9//FEvv/yyZsyYIUl66qmnlJKSopo1a173mszTeeoybae2H80psxxxUcHq1zJW9cP9y2xNAAAAOD4GPqCc2e12TZw4Ua+//rpycnIUEhKiYcOG6Z///OdVbd+s7Vl6Z/ouFZbYyjyHp5tZfZ+zqkVcRJmvDQAAAMfEwAdUkKysLLVt21YLFy6UJLVs2VLDhw9X1apVJUkTNx/RuzN3l2sGk6RPnonR8/fWKtd9AAAA4BgY+IAKZLfbNWbMGHXu3Fl5eXkKDw9XSkqKbFFN1enLHRWSwSRp0N+b0PQBAABUAgx8gAEOHTqkpKQkrVq1Su6hkYpsN0x2U8XdTdPTzaz5HR7id/oAAABcHAMfYBCbzabkz4eo344SedZoWObrn57+sfIPbJQk+VkeU1jzN694Pb5WsKa2/z137wQAAHBht/80aAClYjab5d/06Z+HvTJ+3+Xi/o2Xh73r2XYkR+PWHSrTfQEAAOBYGPgAgxQUlWjYisyfv7mJZ/PdLFthvn5aOkomL1+5V4m84bnDVmSqoKikzPYGAACAY2HgAwwyZ9fx23qo+vXkrB6vknM/KvihF+TmG3zjc/OLNHfXiTLPAAAAAMfAwAcYZPzGw1d8n70yTYf7NNfR5OdVkpd9zWtOTempw32a60RqR9lLiq96/dLJTJ3fOkce4fUUEN/spnJ88ascAAAAcB0MfIABjmZf1M5juVccC37oBXlWu0O2/HM6M3/wVdec2zpHBYe2yeTupbCn35bJzf2K1+12m35aOEyy21Xlz6/IZL65u37uPJajo9kXS//DAAAAwGEx8AEG2P2rYU+STG7uPw9y7l4qOLhV57fNu/xa0dljylmRJkkKfiRRHmFRV11/futcFZ48IH/r4/KKaHRLefZkXZ0HAAAAzo+BDzDA7usMWB5hUQp+JFGSlL18nIrOHpPdVqIzcwbIXnxJ3nXiFNC0+VXXFZ8/o5zVX8jsE6jgRxPLLA8AAACcGwMfYIC9x89d97WAps3lXTde9uJLOjNngHLWjFfhyQMyeweoSrNOMl3jjp7ZS0bLXpiv4EcS5OYTeMt59mRdPw8AAACcFwMfYIDc/MLrvmYymVSlWSeZfQJVePKAzm2YKkkK/fOrcg+octX5FzM36+L+9fKs2VD+1idKmafs7xYKAAAA4zHwAQYoLLnxg9bd/UMV/HCry9/73vWg/Bo9dNV5tsIC/bR4hGQyq8oTr16z/bu5PLZSXQcAAADH5v7bpwAoa55uNx7M7LYSXdiz/PL3hacOylZYILOn9xXnnds0XSXnfpR/7J/lHlpTtsL8Xy/0f38vufyaycP7qsHQ0433fgAAAFwRAx9ggCAfzxu+fm7jNF06tk8mLz+ZvXxVnH1c2cvHqMqTr19xXnHuaUlS3s5Fytu56LrrXdi7Uhf2rpQkRbw8Vu7B1X6Vx6MUPwUAAAAcHW/rAwaIrnn9G6tcOpmpnLWTJEmhj7dXWLM3JZmUt2Oh8r/fUi55LBG3fqMXAAAAOD4aPsAAMRFB1zxuK7qks3MGSLZi+TZ8QP6WP0qSAu75m85vmaWz8z9XjTZD5eb78/Vhzd9UWPM3r7vPyQnddOnoHvlZHrvhedfLAwAAAOdGwwcYICby2gNWzso0FZ09Kjf/UIU++drl4yEPJ8gjrJZKLmTr7MKhZZ7HwsAHAADgkhj4AANEhfgq9ldDX/6h7Tq/da4kqcpTHa94np7J3UNhT78tubkrf/8G5e1aWmZZYiODFRXiW2brAQAAwHEw8AEGeeH+2pe/LinI09n5gyXZ5R/fTD71ml51vme1egp+8HlJ0k/LRqs451SZ5HjxFzkAAADgWkx2u/3GDwQDUC4Kikp0/6fLlGPgQ8+DfTy0sftj8vZwMywDAAAAyg8NH2AQbw83vfZofUP2/u/7PPcF5MjLnf8MAAAAuCr+Tw8wUOsH6iouKrjC9zWZTCo49o1Gv/Uv/e1vf9OJEycqPAMAAADKHwMfYCA3s0n9WsbK061i/1X0dDPr7YeqKSgwQHPmzFF0dLQmTJggPuENAADgWhj4AIPVD/dX3+esMlXQfiZJfZ+z6u2XXtSePXv05JNPKjs7Wy+88IKeffZZnTpVNjeDAQAAgPEY+AAH0CIuQp88E1PuQ59J0ifPxKhFXIQkKTIyUvPnz9eYMWMUEBCgWbNmKTo6WlOmTKHtAwAAcAHcpRNwILO2Z+md6btUWGIr87W93M367Fnr5WHv144cOaK2bdtqyZIlkqSWLVtq+PDhqlq1aplnAQAAQMVg4AMcTObpPHWZtlPbj+aU2ZrxtYLVr2Ws7qjqf8Pz7Ha7UlJS9NZbbykvL09Vq1bV8OHD1bJlyzLLAgAAgIrDwAc4oBKbXePWHdKwFZm39Zy+YB8PvfZofbV+oK7czDf/gdEffvhBbdq00fLlyyVJ//jHPzR06FCFhYWVOgsAAAAqHgMf4MAKiko0d9cJfbHxsHYeu/nGLzYyWK1+V1vNYmqU+qHqNptNI0eOVJcuXXTx4kWFh4dr1KhRatGiRanWAwAAQMVj4AOcxNHsi9qTlavdWbnak3VOuflFKiyxydPNrCAfD1kiAhUTESRLRJCiQnzLbN+DBw+qdevWWrVqlSTp3//+tz7//HOFhoaW2R4AAAAoHwx8AH6TzWbTsGHD1LVrV+Xn56t69eoaPXq0nn76aaOjAQAA4AYY+ADctMzMTCUlJWnt2rWSpFatWmnw4MEKCQkxOBkAAACuhefwAbhp9evX18qVKzVo0CB5e3srIyNDFotF8+fPNzoaAAAAroGGD0CpfPfdd0pKStKGDRskSUlJSRo0aJCCgoIMTgYAAID/ouEDUCoNGzbUmjVr1K9fP3l5eSk1NVUWi0WLFi0yOhoAAAD+Dw0fgNv2zTffKDExUZs3b5YktWvXTv3791dgYKDByQAAACo3Gj4At61Ro0Zat26d+vTpI09PT6WkpCgmJkZLly41OhoAAEClRsMHoEzt3btXCQkJ2rp1qyTp5ZdfVt++fRUQEGBwMgAAgMqHhg9AmYqOjtaGDRv08ccfy8PDQyNHjpTVatWKFSuMjgYAAFDp0PABKDe7du1SYmKitm/fLkl6/fXX1adPH/n5+RmcDAAAoHKg4QNQbqxWqzZt2qQPP/xQ7u7uGjp0qKxWq1avXm10NAAAgEqBhg9AhdixY4cSEhK0a9cumUwmdejQQb1795avr6/R0QAAAFwWDR+ACtGkSRNt2bJFPXr0kNlsVnJyspo0aaJ169YZHQ0AAMBl0fABqHBbt25VYmKi9uzZI5PJpM6dO+ujjz6Sj4+P0dEAAABcCg0fgArXtGlTff3113r33XdlMpk0YMAAxcXFaePGjUZHA9pi9rgAACAASURBVAAAcCk0fAAMtWXLFiUkJOibb76R2WzW22+/rQ8//FDe3t5GRwMAAHB6NHwADHXPPfdo27Zt6tq1qySpb9++io+P15YtWwxOBgAA4Pxo+AA4jI0bNyoxMVHfffedzGazunbtql69esnLy8voaAAAAE6JgQ+AQ8nPz1ePHj00cOBA2e12WSwWpaWlqWnTpkZHAwAAcDoMfAAc0rp165SYmKjMzEy5ubnp3Xff1fvvvy9PT0+jowEAADgNBj4ADuvixYt67733lJycLLvdLqvVqvT0dDVp0sToaAAAAE6BgQ+Aw1u9erWSkpJ08OBBubu76/3339e7774rDw8Po6MBAAA4NAY+AE7hwoUL6tatm4YOHSpJiouLU1pamqxWq8HJAAAAHBePZQDgFPz8/DRkyBAtX75cderU0fbt23X33Xfrk08+UXFxsdHxAAAAHBINHwCnk5eXp3feeUcjRoyQJDVt2lTp6emKjo42OBkAAIBjoeED4HT8/f01fPhwLVmyRLVq1dLWrVsVHx+vPn360PYBAAD8Ag0fAKd27tw5denSRaNHj5Yk3XvvvUpLS1OjRo0MTgYAAGA8Gj4ATi0wMFCjRo3SwoULFRkZqc2bNysuLk79+vVTSUmJ0fEAAAAMRcMHwGXk5uaqc+fOGjdunCTpd7/7nVJTU9WwYUODkwEAABiDhg+AywgKCtLYsWM1b9481axZUxs2bFCTJk00aNAg2j4AAFAp0fABcEnZ2dl68803lZ6eLkl68MEHlZqaqvr16xucDAAAoOLQ8AFwSSEhIUpLS9NXX32l6tWra+3atbJarRoyZIhsNpvR8QAAACoEDR8Al/fTTz+pQ4cOmjBhgiTp4Ycf1rhx41SvXj2DkwEAAJQvGj4ALi80NFTjx4/XzJkzFR4erlWrVslqtWr48OG0fQAAwKXR8AGoVM6cOaM33nhDkydPliT98Y9/1NixY1WnTh1jgwEAAJQDGj4AlUpYWJgmTZqkqVOnKiwsTMuXL1dMTIxGjRol3v8CAACuhoYPQKV1+vRpvfbaa5o2bZok6fHHH9eYMWNUq1Ytg5MBAACUDRo+AJVWeHi4pk6dqilTpqhKlSpasmSJLBaLxo4dS9sHAABcAg0fAEg6deqUXnnlFc2cOVOS9OSTTyolJUWRkZEGJwMAACg9Gj4AkFStWjVNnz5dEyZMUEhIiBYuXCiLxaK0tDTaPgAA4LRo+ADgV06cOKH27dtrzpw5kqRmzZpp9OjRqlmzpsHJAAAAbg0NHwD8So0aNTR79mxlZGQoODhY8+bNU3R0tL744gvaPgAA4FRo+ADgBrKysvTSSy9p/vz5kqS//vWvGjVqlKpXr25wMgAAgN9GwwcANxAREaG5c+dq3LhxCgwM1FdffaXo6GhNmjSJtg8AADg8Gj4AuElHjx5Vu3bttGjRIknSs88+qxEjRig8PNzgZAAAANdGwwcANykqKkoLFixQSkqKAgICNGPGDEVHR+vLL780OhoAAMA10fABQCkcPnxYbdq00bJlyyRJ/+///T8NGzZMVatWNTgZAADA/9DwAUAp1K5dW0uWLNGIESPk5+enqVOnKjo6WtOnTzc6GgAAwGU0fABwmw4dOqQ2bdpoxYoVkqR//etfGjJkiKpUqWJwMgAAUNnR8AHAbapbt66WLl2qoUOHytfXV5MmTVJ0dLRmz55tdDQAAFDJ0fABQBn6/vvv1bp1a61evVqS9MILLyg5OVmhoaEGJwMAAJURDR8AlKE77rhDK1asUHJysnx8fDR+/HhZLBbNnTvX6GgAAKASouEDgHJy4MABJSUlad26dZKkhIQEDR48WMHBwQYnAwAAlQUNHwCUkzvvvFOrVq3SgAED5O3trfT0dFksFi1YsMDoaAAAoJKg4QOACvDdd98pMTFRGzdulCS1adNGAwYMUFBQkMHJAACAK6PhA4AK0LBhQ61du1Z9+/aVl5eXxo4dq5iYGC1ZssToaAAAwIXR8AFABdu3b58SExO1ZcsWSdJLL72k/v37KyAgwOBkAADA1dDwAUAFa9y4sdavX69PP/1Unp6eGj16tGJiYrRs2TKjowEAABdDwwcABtqzZ48SExO1detWSdKrr76qzz77TP7+/gYnAwAAroCGDwAMZLFYtGHDBn300Ufy8PDQ8OHDZbVatXLlSqOjAQAAF0DDBwAOYteuXUpISNCOHTskSW+88YY+/fRT+fn5GZwMAAA4Kxo+AHAQVqtVmzdv1gcffCB3d3cNGTJEsbGxWrNmjdHRAACAk6LhAwAHtH37diUkJGj37t0ymUzq1KmTPv74Y/n6+hodDQAAOBEaPgBwQHFxcfr666/1/vvvy2w2a9CgQWrSpInWr19vdDQAAOBEaPgAwMF9/fXXSkxM1N69e2UymfTWW2/pP//5j3x8fIyOBgAAHBwNHwA4uLvvvltbt25V9+7dZTKZ1L9/f8XHx2vTpk1GRwMAAA6Ohg8AnMjmzZuVkJCgb7/9VmazWV26dNEHH3wgb29vo6MBAAAHRMMHAE7k3nvv1fbt29WlSxdJ0meffaamTZtqy5YtBicDAACOiIYPAJzUhg0blJiYqP3798vNzU3dunVTjx495OXlZXQ0AADgIBj4AMCJ5efn6/3339egQYNkt9sVExOjtLQ0xcfHGx0NAAA4AAY+AHABa9euVVJSkjIzM+Xu7q53331X7733njw9PY2OBgAADMTABwAu4uLFi3r33XeVnJwsSYqNjVV6erpiY2MNTgYAAIzCTVsAwEX4+vpq8ODBWrlyperVq6edO3fq7rvv1kcffaSioiKj4wEAAAPQ8AGAC8rLy1O3bt00bNgwSVJ8fLzS09NlsVgMTgYAACoSDR8AuCB/f38NHTpUy5cvV506dbRt2zbFx8erd+/eKi4uNjoeAACoIDR8AODizp8/r3feeUcjR46UJN19991KT09X48aNDU4GAADKGw0fALi4gIAAjRgxQosXL1ZUVJS+/vprxcXF6bPPPqPtAwDAxdHwAUAlcu7cOb311lsaM2aMJOm+++5TWlqa7rrrLoOTAQCA8kDDBwCVSGBgoFJSUrRw4UJFRERo06ZNatKkiQYMGKCSkhKj4wEAgDJGwwcAlVROTo46d+6s1NRUSdLvf/97paamqkGDBgYnAwAAZYWGDwAqqeDgYI0bN05z585VjRo1tH79ejVp0kTJycmy2WxGxwMAAGWAhg8AoOzsbHXq1EkZGRmSpIceekjjxo1T/fr1DU4GAABuBw0fAEAhISFKT0/X7NmzVb16da1Zs0axsbEaOnQobR8AAE6Mhg8AcIWzZ8+qQ4cOmjhxoiTpkUce0bhx41S3bl2DkwEAgFtFwwcAuEKVKlU0YcIEzZgxQ+Hh4Vq5cqViYmI0YsQI2j4AAJwMDR8A4LrOnDmj119/XVOmTJEkPfbYYxo7dqxq165tcDIAAHAzaPgAANcVFhamyZMn68svv1RYWJiWLVsmi8WilJQU8X4hAACOj4YPAHBTTp8+rVdffVXTp0+XJD3xxBMaM2aMoqKiDE4GAACuh4YPAHBTwsPDNXXqVE2aNEmhoaFavHixLBaLxo0bR9sHAICDouEDANyykydP6uWXX9bs2bMlSX/5y1+UkpKiiIgIg5MBAIBfouEDANyy6tWra+bMmRo/frxCQkK0YMECRUdHKz09nbYPAAAHQsMHALgtx48fV/v27TV37lxJUvPmzTV69GjVqFHD4GQAAICGDwBwW2rWrKmvvvpK6enpCgoK0ty5cxUdHa0JEybQ9gEAYDAaPgBAmcnKylK7du20YMECSVKLFi00cuRIVatWzeBkAABUTjR8AIAyExERoXnz5mns2LEKDAzUrFmzFB0drcmTJ9P2AQBgABo+AEC5OHr0qNq2bavFixdLkp577jkNHz5c4eHhBicDAKDyoOEDAJSLqKgoLVy4UKNHj5a/v7+mT5+u6OhoTZ061ehoAABUGjR8AIByd/jwYbVp00bLli2TJP3jH//Q0KFDFRYWZnAyAABcGw0fAKDc1a5dW0uWLNGIESPk5+enKVOmKDo6WjNnzjQ6GgAALo2GDwBQoQ4dOqTWrVtr5cqVkqTnn39en3/+uapUqVIm6x/Nvqjdx3K1OytXe4+fU25+oQpL7PJ0MynIx1PRNQMVExGkmMggRYX4lsmeAAA4KgY+AECFs9lsGj58uLp27aqLFy+qevXqGjVqlP7617+War2CohLN2XVc4zce1s5juTd9XWxkkF68v46aW2vI28OtVHsDAODIGPgAAIb5/vvvlZSUpDVr1kiSXnzxRSUnJyskJOSmri+x2TVu3SENW5GpnPyiUucI9vHQa4/WV+sH6srNbCr1OgAAOBoGPgCAoWw2m4YMGaLu3bsrPz9fNWvW1OjRo9WsWbMbXpd5Ok9dpu3U9qM5ZZYlLipY/VrGqn64f5mtCQCAkRj4AAAOYf/+/UpKStL69eslSUlJSRo4cKCCg4OvOnfW9iy9M32XCktsZZ7D082svs9Z1SIuoszXBgCgonGXTgCAQ2jQoIFWr16t/v37y8vLS6mpqYqJidGiRYuuOG/i5iPq9OWOchn2JKmwxKY3v9yhiZuPlMv6AABUJBo+AIDD+fbbb5WYmKhNmzZJktq2basBAwZo+ffn1enLHRWSwSRp0N+b0PQBAJwaAx8AwCEVFxdr4MCB6tGjhwoLCxUVfbc8//qBiivwTy1PN7Pmd3iI3+kDADgtBj4AgEPbt2+fWiUkKqvRP+QVcVep1yk+f0YX9q1S4fH9KjpzRCUXc2W7dEFmLz95hNWSb8MHFNDkSZncPa64Lr5WsKa2/z137wQAOCUGPgCAwxu5KlN9Fn4nu90uk6l0g9eFfat05qt+/zvg5iGTu4fsly5ePuRRJUrh//xI7gFhV1z73lON1O6heqXaFwAAIzHwAQAcWkFRie7/dNltPWdPkgqOfaOCQ9vkXStGHuF15eYTIEmyXbqoC9+sVvbysbIX5su7tlXV/tX7imuDfTy0sftjPJwdAOB0uEsnAMChzdl1/LaHPUnyjmyk4If+Le/a1svDniSZvXwV0ORJhTzWTpJUcHiXis/9eMW1OflFmrvrxG1nAACgojHwAQAc2viNh6/4Pntlmg73aa6jyc+rJC/7mtecmtJTh/s014nUjrKXFN/UPl417rz8dUneT1e9/sWvcgAA4AwY+AAADuto9kXtPJZ7xbHgh16QZ7U7ZMs/pzPzB191zbmtc1RwaJtM7l4Ke/ptmdzcb2qvS1nfXP7aPSj8qtd3HsvR0eyLVx0HAMCRMfABABzW7l8Ne5JkcnP/eZBz91LBwa06v23e5deKzh5Tzoo0SVLwI4nyCIu64fr2kiIV5ZzUuc0zlb18rCTJ964H5eYXcs3z92RdnQcAAEd2c297AgBggN3XGbA8wqIU/EiispeOUvbycfKuHSv3kBo6M2eA7MWX5F0nTgFNm1933ayR7VSc8+vfyTPJ964HVeWpjjfM8xdLjdL8KAAAGIKGDwDgsPYeP3fd1wKaNpd33XjZiy/pzJwBylkzXoUnD8jsHaAqzTrd8PENZt9Amf2CZfL0uXzMt9GDCn64lcy/OPZre7KunwcAAEdEwwcAcFi5+YXXfc1kMqlKs046MfZ1FZ48oMKTByRJoX9+Ve4BVW64bo1WAy5/XXIhW3k7Fyt3w1TlH9ikKs07y++uB6+T5/bvFgoAQEWi4QMAOKzCkhs/KtbdP1TBD7e6/L3vXQ/Kr9FDt7SHm1+Ign7/D4X9tYvsxYU6O2+Qis+fuU4e2y2tDQCA0Rj4AAAOy9Pt+h/LlCS7rUQX9iy//H3hqYOyFRaUai/fO++TW2C47EWXdHHf6uvk4Y9NAIBz4U8uAIDDCvLxvOHr5zZO06Vj+2Ty8pNbYFUVZx9X9vIxpd7vvx8FLco5eZ08HqVeGwAAIzDwAQAcVnTNwOu+dulkpnLWTpIkhT7eXmHN3pRkUt6Ohcr/fkup9ivOPSVJMnt6X/N1S8T18wAA4IgY+AAADismIuiax21Fl3R2zgDJVizfhg/I3/JHede2KuCev0mSzs7/XCUXr3ykg91WcsO98vauUEneT5Ikr8joW8oDAICjYuADADismMhrD1g5K9NUdPao3PxDFfrka5ePhzycII+wWiq5kK2zC4decc2pCV2Vu2GqCs8cuWL4K849rZy1E3V2frIkybN6ffnUv+ea+1oY+AAATobHMgAAHFZUiK9iI4O089j/2rr8Q9t1futcSVKVpzrKzed/H7M0uXso7Om3dSKjs/L3b1DerqXyt/5JklR8/iflrEpXzqp0yewus5ev7MWFshf97yYvnjXuVPhzPWUyXf1+aGxksKJCfMvrRwUAoFzQ8AEAHNoL99e+/HVJQZ7Ozh8syS7/+Gbyqdf0qvM9q9VT8IPPS5J+WjZaxTk//15eWLNOCry/pbwi7pKbX7BshfmS3S63oGrybfB7hf21i6q/2F9u/iHXzPHiL3IAAOAsTHa7/cYPOQIAwEAFRSW6/9NlyjHwoefBPh7a2P0xeXu4GZYBAIDSoOEDADg0bw83vfZofUP2/u97og9VuciwBwBwSgx8AACH1/qBuoqLCq7wfU0mkwqOfaOhbzynv//97zp9+nSFZwAA4HYw8AEAHJ6b2aR+LWPl6Vaxf2x5upn12t2B8vP10dSpUxUdHa3JkyeL34YAADgLBj4AgFOoH+6vvs9ZZaqg/UyS+j5nVY+OL2nPnj3605/+pDNnzuhf//qXnnnmGZ04caKCkgAAUHoMfAAAp9EiLkKfPBNT7kOfSdInz8SoRVyEJKlOnTpavHixUlJSFBgYqNmzZys6OloZGRm0fQAAh8ZdOgEATmfW9iy9M32XCktsZb62l7tZnz1rvTzs/drRo0fVvn17LViwQJL01FNPadSoUYqMjCzzLAAA3C4GPgCAU8o8nacu03Zq+9GcMlszvlaw+rWM1R1V/W94nt1uV0ZGhjp16qScnBwFBgZq4MCBat26tUymivrQKQAAv42BDwDgtEpsdo1bd0jDVmTe1nP6gn089Nqj9dX6gbpyM9/8wHb8+HG98sor+uqrryRJjz/+uFJSUlS7Ng9pBwA4BgY+AIDTKygq0dxdJ/TFxsPaeezmG7/YyGC1+l1tNYupUern7Nntdk2ePFlvvPGGzp49K39/f/Xt21ft27eX2cyvygMAjMXABwBwKUezL2pPVq52Z+VqT9Y55eYXqbDEJk83s4J8PGSJCFRMRJAsEUGKCvEts31PnTql119/XdOmTZMkPfLIIxo7dqzq1atXZnsAAHCrGPgAAChD06ZN06uvvqoff/xRvr6+6t27t9544w3aPgCAIRj4AAAoY2fOnFGHDh00adIkSdIDDzygcePGqUGDBgYnAwBUNrzdCABAGQsLC9PEiRM1a9Ys1ahRQ+vWrVNsbKz69++vkpISo+MBACoRGj4AAMpRdna23nzzTaWnp0uS7r33XqWmpqpx48YGJwMAVAY0fAAAlKOQkBClpaVp/vz5ioyM1ObNmxUXF6fevXurqKj0j5IAAOBm0PABAFBBcnNz1aVLF6WkpEiS4uPjlZqaKqvVanAyAICrouEDAKCCBAUFafTo0VqyZIlq166tbdu2qWnTpvrggw9UWFhodDwAgAui4QMAwADnz59X9+7dNWzYMEmS1WpVamqq4uPjDU4GAHAlNHwAABggICBAQ4cO1cqVK3XHHXdo165duvfee/Xee+/p0qVLRscDALgIGj4AAAx24cIFvf/++0pOTpbdblejRo2Umpqq++67z+hoAAAnx8AHAICDWLdunVq3bq39+/fLbDarc+fO+s9//iMfHx+jowEAnBQf6QQAwEE88MAD2rFjh9555x1JUv/+/dWkSROtW7fO4GQAAGdFwwcAgAPavHmzkpKStG/fPplMJnXo0EGffPKJ/Pz8jI4GAHAiNHwAADige++9V9u2bdN7770ns9ms5ORkWa1WrVy50uhoAAAnQsMHAICD27Ztm5KSkrRr1y5J0iuvvKLPPvtMAQEBBicDADg6Gj4AABxcfHy8tmzZog8//FAeHh4aMWKEYmJitGTJEqOjAQAcHA0fAABOZPfu3UpKStLWrVslSW3btlX//v0VFBRkcDIAgCOi4QMAwInExMRo48aN6t27tzw9PTVmzBhFR0dr/vz5RkcDADggGj4AAJzUvn371Lp1a23atEmS1KpVKw0aNEihoaEGJwMAOAoaPgAAnFTjxo21bt069e/fX97e3srIyFB0dLRmz55tdDQAgIOg4QMAwAXs379fbdq00dq1ayVJ//znPzVkyBCFhYUZnAwAYCQaPgAAXECDBg20atUqff755/L19dXkyZPVuHFjTZ061ehoAAAD0fABAOBiDh48qLZt22rFihWSpOeee07Dhg1TtWrVDE4GAKhoNHwAALiYevXqaenSpRoxYoT8/f01ffp0NW7cWBMmTBDv8wJA5ULDBwCACzty5IjatWunxYsXS5KefvppjRw5UjVr1jQ4GQCgItDwAQDgwmrVqqWFCxdq7NixCgoK0pw5c9S4cWOlpqbS9gFAJUDDBwBAJZGVlaX27dtr3rx5kqQnn3xSo0ePVlRUlMHJAADlhYYPAIBKIiIiQnPmzFFGRoZCQkK0cOFCRUdHa/To0bR9AOCi/n97dx4dVX3/f/w12UjINiRhDWEzskw2EnbUqq0//bZYqy31e2rVSgJqga8sFRSwlFaUTSBACASSCcgmEAQBkS0skR2BEJaioojsKiZBIOvM/P7AUpXVkOTOTJ6PczgnM3fu577yB+fw5nXnfmj4AACogc6ePavevXtr6dKlkqRf/epXSk9PV7NmzYwNBgCoVDR8AADUQA0aNNCSJUv0zjvvKCwsTNnZ2YqOjtbUqVNlt9uNjgcAqCQ0fAAA1HBff/21+vbtq0WLFkmSfvGLXygjI0ORkZEGJwMA3CkaPgAAari6detq4cKFWrJkierXr6+cnBzFxsZq4sSJstlsRscDANwBGj4AAHDV+fPn1b9/f82dO1eS1KVLF1mtVrVu3drgZACAiqDhAwAAV4WGhmrOnDlasWKFGjVqpO3bt6tt27YaM2aMysvLjY4HAPiZaPgAAMB1FRQUaODAgcrMzJQkdejQQVarVdHR0QYnAwDcLho+AABwXWazWVarVatXr1ZERIR2796thIQEjRw5UmVlZUbHAwDcBho+AABwSxcuXNDgwYOVlpYmSWrbtq0yMzPVtm1bg5MBAG6Ghg8AANxSUFCQpk+fruzsbDVr1ky5ubnq0KGDhg8frtLSUqPjAQBugIYPAAD8LBcvXtTQoUM1ZcoUSVJ0dLQyMzPVvn17g5MBAH6Khg8AAPwsAQEBmjx5snJychQZGamDBw+qU6dOevXVV1VcXGx0PADAD9DwAQCACrt8+bKGDx+uCRMmyOFwqHXr1rJarerSpYvR0QAAYuADAACVYMeOHerRo4eOHDkik8mkAQMG6PXXX1ft2rWNjgYANRq3dAIAgDvWuXNn7du3T6+++qpMJpMmTJiguLg45eTkGB0NAGo0Gj4AAFCpPvroI/Xo0UMHDx6UJPXt21ejRo1SQECAwckAoOah4QMAAJWqffv22rNnj4YPHy4vLy+lpKQoNjZWGzZsMDoaANQ4NHwAAKDK5ObmqkePHsrNzZUkvfDCCxo7dqyCgoIMTgYANQMNHwAAqDJt27bVrl279Prrr8vb21tpaWmKjo7WmjVrjI4GADUCDR8AAKgWBw8eVGJionbv3i1J6tGjhyZMmCCz2WxwMgBwXzR8AACgWkRHR2vbtm0aM2aMatWqpczMTEVFRWnlypVGRwMAt0XDBwAAqt2RI0eUmJio7du3S5KefvppJScnKzQ01OBkAOBeaPgAAEC1a926tT788ENNmDBBfn5+mjt3rqKiovTuu+8aHQ0A3AoNHwAAMNTRo0eVlJR0dZP2J598UikpKapbt67ByQDA9dHwAQAAQ0VGRmrjxo1KSUmRv7+/Fi1aJIvFooULF4r/lwaAO0PDBwAAnMaxY8fUq1cvZWdnS5KeeOIJpaamqkGDBgYnAwDXRMMHAACcRvPmzbVu3TqlpaUpMDBQS5culcVi0Zw5c2j7AKACaPgAAIBTOnHihJ5//nmtXr1aktStWzelpaUpPDzc4GQA4Dpo+AAAgFOKiIjQqlWrlJmZKbPZrPfff18Wi0UZGRm0fQBwm2j4AACA0zt9+rRefPFFrVixQpL08MMPa8aMGWratKnByQDAudHwAQAAp9eoUSO99957mjdvnkJCQrR27VpFR0dr2rRpstvtRscDAKdFwwcAAFzKuXPn1KdPHy1ZskSS9OCDDyo9PV0tWrQwOBkAOB8aPgAA4FLq16+vrKwsLV68WHXr1tXGjRsVExOjyZMn0/YBwE/Q8AEAAJf1zTff6KWXXtKCBQskSffee68yMjLUsmVLg5MBgHOg4QMAAC4rLCxM8+fP17Jly9SgQQNt2bJFcXFxeuutt2Sz2YyOBwCGo+EDAABu4dtvv9WAAQP09ttvS5I6deokq9Uqi8VicDIAMA4NHwAAcAshISGaPXu23n//fYWHh2vnzp2Kj4/XqFGjVF5ebnQ8ADAEDR8AAHA7hYWFevnll5Weni5JateunaxWq2JjYw1OBgDVi4YPAAC4neDgYM2cOVNr165VkyZNtGfPHrVv317//Oc/VVpaanQ8AKg2NHwAAMCtfffdd3rllVc0bdo0SVJsbKwyMzOVkJBgcDIAqHo0fAAAwK0FBgYqNTVVGzduVIsWLZSXl6eOHTvqtddeU0lJidHxAKBK0fABAIAa49KlSxo2bJgmT54sh8Mhi8WizMxMdezY0ehoAFAlaPgAAECN4e/vr+TkZH344Ydq2bKlDh8+rC5dumjw4MEqKioyOAylAAAAIABJREFUOh4AVDoaPgAAUCMVFRVp+PDhmjBhgux2u1q2bKnMzEx17drV6GgAUGkY+AAAQI22c+dOJSYm6vDhwzKZTOrXr59Gjhwpf39/o6MBwB3jlk4AAFCjderUSXv37tXQoUPl4eGh5ORkxcXFadOmTUZHA4A7RsMHAADwvb1796pHjx7Ky8uTJPXu3VujR49WYGCgwckAoGJo+AAAAL6XkJCg3bt3a8SIEfLy8lJqaqpiYmK0fv16o6MBQIXQ8AEAAFxHXl6eevToob1790qSevbsqbfeekvBwcEGJwOA20fDBwAAcB2xsbHauXOn3nzzTfn4+Cg9PV3R0dH64IMPjI4GALeNhg8AAOAWDh8+rMTERO3cuVOS9Oyzzyo5OVl16tQxOBkA3BwNHwAAwC1YLBZt3bpV48aNk6+vr95++21ZLBa99957RkcDgJui4QMAAPgZPvnkEyUmJmrr1q2SpD/96U+aPHmywsLCDE4GANei4QMAAPgZWrZsqZycHE2aNEm1a9fWggULFBUVpaysLKOjAcA1aPgAAAAq6LPPPlPPnj2vbtLevXt3paSkqH79+sYGA4Dv0fABAABU0F133aXs7GylpqYqICBAWVlZioqK0vz588X/qQNwBjR8AAAAleD48ePq1auX1q1bJ0l67LHHNH36dDVs2NDgZABqMho+AACAStC0aVOtWbNG6enpCgoK0vLly2WxWDR79mzaPgCGoeEDAACoZCdPntQLL7ygVatWSZJ+/etfKy0tTREREQYnA1DT0PABAABUssaNG2vlypWaPXu2zGazPvjgA0VFRWnmzJm0fQCqFQ0fAABAFTpz5ox69+6tZcuWSZIeeughzZw5U82aNTM2GIAagYYPAACgCjVs2FDvvvuuFixYoNDQUK1fv17R0dGaOnWq7Ha70fEAuDkaPgAAgGry1VdfqW/fvlq8eLEk6f7771d6eroiIyMNTgbAXdHwAQAAVJN69epp0aJFysrKUr169bR582bFxsYqOTlZNpvN6HgA3BANHwAAgAHOnz+vfv36ad68eZKkrl27ymq1qlWrVgYnA+BOaPgAAAAMEBoaqrlz52r58uVq1KiRtm3bpri4OI0dO1bl5eVGxwPgJmj4AAAADFZQUKCBAwcqMzNTktShQwdlZmYqKirK4GQAXB0NHwAAgMHMZrOsVqs++OADNW7cWLt371Z8fLxGjhypsrIyo+MBcGE0fAAAAE7kwoULGjx4sNLS0iRJ8fHxyszMVFxcnMHJALgiBj4AAAAnlJ2drZ49e+qLL76Ql5eXhg4dqmHDhsnHx6dS1j+Rf1kHThbqwKlCHTp9QYVFpSq1OeTjaVKwn4+iGgUpJjxYMY2DFVGndqVcE0D1Y+ADAABwUhcvXtSQIUOUkpIiSYqJiVFmZqbatWtXofWKy2xakXdac3cc1/6Thbd9XlzjYD3TuZkejW0oX2/PCl0bgDEY+AAAAJxcTk6OEhMT9dlnn8nT01ODBg3SP/7xD/n6+t7W+Ta7Q9atxzR141EVFFX8O4FmP2/1eTBSifc0l6eHqcLrAKg+DHwAAAAu4PLly/r73/+uiRMnyuFwqE2bNrJarercufNNzzv61UUNytqvfScKKi1LfIRZ47rHKbJeQKWtCaBqMPABAAC4kO3btysxMVFHjhyRh4eHBgwYoNdff11+fn7XfHbZvlMavCRPpTZ7pefw8fTQ2D/E6vH48EpfG0DlYeADAABwMcXFxRoxYoTGjRsnu92uu+++WxkZGbrvvvuufmb+ri81dOmBKs1hkvTGEzF6qmOTKr0OgIpj4AMAAHBRu3fvVo8ePXTo0CGZTCb17dtXo0aN0rpPCtR/UW61ZDBJmvhkW5o+wEkx8AEAALiwkpISvfHGGxo1apTKy8vVLLazPB99TeWVfxfnDfl4emjVS/fxnT7ACTHwAQAAuIHc3Fw91yNR52L+rFrhrSu8jr3ksi5/ukPFx/ap5Mynsl34WpJDngEhqhURraB2v5VPg7uuOS+hiVmLX+jK0zsBJ8PABwAA4Camb/pUo9d8IofDIZOpYoPXqbReKs8/c/W1ybuW5JAc5SXfv+GhOg88p6BOv7/m3GG/aaNe97Wo0HUBVA0vowMAAADgzhWX2TQ955gkVXjYkySHzSaf+ncpIO5h+d3VQV7B9eRw2FX29XF9u36mSr7MU/5Gq7zDIuR3V4cfnTt141E907kpm7MDTsTD6AAAAAC4cyvyTt/Rpur/EfbbgWrYY5ICE7rJK7ieJMlk8pBPveaq98d/yDs0QpJUuPPda84tKCrTyrwz17wPwDgMfAAAAG5g7o7jP3qdv2mWjo9+VCcmPSXbxfzrnnNu4XAdH/2ozmT2k8NWLknyjYi+4TU8vGupdpsrWz+Unj163c/M+UkOAMZi4AMAAHBxJ/Iva//Jwh+9Z77vafnUv0v2ogv6ZlXyNedc2LNCxcf2yuRVS2G/fVkmz9v7po+HX+CVHxzXfwzo/pMFOpF/+ef9AgCqDAMfAACAizvwk2FPkkyeXlcGOa9aKv58j77b+/7VY2XnT6pg4yxJkvmB5+QdFnHb1yr58qAkyTus6Q0/c/DUtXkAGIOBDwAAwMUduMGA5R0WIfMDz0mS8jdYVXb+pBx2m75ZMV6O8hL5NotXYLtHb/s6JWeP6vIn2yVJAbEP/ew8AKofAx8AAICLO3T6wg2PBbZ7VL7NE+QoL9E3K8ar4MO5Kj37qTx8AxXarf9tP9HTXnJZ36x4S3LYv3+K5yM3/OzBUzfOA6B6MfABAAC4uMKi0hseM5lMCu3WXx5+QSo9+6kubF8sSQp5pLe8AkNva32H3aZvlo9T+fmT8qjlr7DfDZbJ48ZbLxRWwtNCAVQOBj4AAAAXV2pz3PS4V0CIzPc/e/V17db3yv/7p23eisNh1/n3J6ros90yeddS3e7D5R0Sfos813+gC4Dqx8AHAADg4nw8b35bpsNu06WDG66+Lj33ueylxbdc1+Fw6Ns1qbp0aJPk6aW6vx8m34io28jDPzEBZ8HfRgAAABcX7Odz0+MXdmSp5ORhmWr5yzOorsrzTyt/Q/ot181fP0MXc1dLHp6q+9hg+TVPuM083rf1OQBVj4EPAADAxUU1CrrhsZKzR1WwZYEkKeT/vaCwbgMkmXQxd7WKPtt9w/PyN83Sd3tWSCYPhXUboNqtut52nujwG+cBUL0Y+AAAAFxcTHjwdd+3l5Xo/Irxkr1ctVvdo4DoX8q3aawCO/xOknR+1WTZLl+7hULB1gW6sCNLkkkh/9NX/lEPVEoeANWPgQ8AAMDFxTS+/oBVsGmWys6fkGdAiEL+p8/V9+vc/xd5hzWR7VK+zq9O+dE5F3a/p8IP50m60ggGxj38s/NEM/ABTsPL6AAAAAC4MxF1aiuucbD2n/xvW1d0bJ++27NSkhT6m37y9PvvbZYmL2+F/fZlnXl7oIo+2a6LeeuvbqSen/39d/tMHirY9o4Ktr1zw+s2/MtEeQXV/dF7cY3NiqhTu7J+NQB3iIYPAADADTzduenVn23FF3V+VbIkhwISusmvRbtrPu9Tv4XM9z4lSfo2e4bKC859f+T7LR4cdtkvFdz0j+zXbr/wzA9yADCeyeFw3HzjFgAAADi94jKbOo/KVoGBm56b/by1Y8iv5Ot9403ZAVQvGj4AAAA34OvtqT4PRhpy7f/0Bw3z81RWfNmQDACuj4EPAADATSTe01zxEeZqv67JZFLJqSNaPekVRUdHa+3atdWeAcD1MfABAAC4CU8Pk8Z1j5OPZ/X+E8/H00Ppve5Xu4R4ffnll3rkkUfUo0cP5efnV2sOANdi4AMAAHAjkfUCNPYPsTJV0/VMksb+IVbd7muvHTt2aPTo0apVq5ZmzZoli8WipUuXVlMSANfDQ1sAAADc0PxdX2rY0gOqyn/omSS98USMnurY5Efvf/zxx0pKStLWrVslSX/84x81ZcoU1a9fvwrTALgeGj4AAAA39FTHJpr4ZNsqu72zlpeHJj7Z9pphT5JatWqlnJwcTZkyRf7+/lq8eLEsFovmzp0rugagetHwAQAAuLGjX13UoKz92neioNLWTGhi1rjucbqrbsAtP/vFF1/o+eef17p16yRJ3bp107Rp0xQREVFpeQDcGAMfAACAm7PZHbJuPaapG4/e0T59Zj9v9XkwUon3NJenx+1/S9DhcGjWrFkaOHCgCgoKFBgYqHHjxqlXr17y8OCGM6AqMfABAADUEMVlNq3MO6M5O45r/8nbb/ziGpv1bJem6hbT8I42VT99+rT69OmjZcuWSZIeeOABpaen66677qrwmgBujoEPAACgBjqRf1kHTxXqwKlCHTx1QYVFZSq12eXj6aFgP29FhwcpJjxY0eHBiqhTu9Ku63A4tHjxYvXt21dff/21/Pz8NHLkSPXr10+enhUfJgFcHwMfAAAAqt0333yj/v37a968eZKkTp06KSMjQ1FRUQYnA9wLN00DAACg2oWFhWnu3LlasWKFwsPDtXPnTiUkJGjkyJEqK6v49wwB/BgNHwAAAAxVWFiowYMHa8aMGZKkuLg4Wa1WJSQkGJwMcH00fAAAADBUcHCw0tLSlJ2drebNm2v//v3q2LGjhgwZouLiYqPjAS6Nhg8AAABO49KlS3rttdc0adIkORwOtWrVShkZGbrnnnuMjga4JAY+AAAAOJ3t27crKSlJ//73v2UymdS3b1+9+eabCgi49WbvAP6LWzoBAADgdLp06aJ9+/Zp2LBh8vDw0JQpUxQTE6N169YZHQ1wKTR8AAAAcGq5ublKTEzUvn37JEmJiYkaP368zGazwckA50fDBwAAAKfWtm1b7dy5U2+++aZq1aolq9Uqi8Wi9957z+hogNOj4QMAAIDLOHLkiJKSkrRt2zZJ0v/+7/9qypQpqlu3rsHJAOdEwwcAAACX0bp1a+Xk5GjSpEmqXbu2Fi5cqDZt2mj+/PmixwCuRcMHAAAAl3Ts2DH16tVL2dnZkqRHH31U06dPV3h4uMHJAOdBwwcAAACX1Lx5c61bt07p6ekKDg7WypUrZbFYNHPmTNo+4Hs0fAAAAHB5p06dUu/evbV8+XJJ0i9/+UvNnDlTLVq0MDgZYCwaPgAAALi88PBwLVu2TAsWLFBYWJg2bNigmJgYJScny2azGR0PMAwNHwAAANzK119/rX79+mnBggWSrmzinpGRoTZt2hicDKh+NHwAAABwK3Xr1tX8+fP13nvvqVGjRtq+fbvatm2rN998U2VlZUbHA6oVDR8AAADcVkFBgQYNGqT09HRJVzZxt1qtio+PNzgZUD1o+AAAAOC2zGazZs6cqXXr1qlZs2bKzc1Vhw4dNGzYMBUXFxsdD6hyNHwAAACoES5evKjXXntNkydPlsPhUOvWrWW1WtWlSxejowFVhoEPAAAANcq2bduUmJiojz/+WCaTSS+99JLeeOMN+fv7Gx0NqHTc0gkAAIAapWvXrsrNzdWQIUPk4eGhSZMmKSYmRtnZ2UZHAyodDR8AAABqrL179yopKUm5ubmSpJ49e2rcuHEym80GJwMqBw0fAAAAaqyEhATt2rVLI0eOlI+Pj9LT0xUVFaUVK1YYHQ2oFDR8AAAAgKTDhw8rKSlJO3bskCT96U9/0qRJk1S3bl2DkwEVR8MHAAAASLJYLNqyZYsmTpyo2rVra8GCBbJYLHrnnXdERwJXRcMHAAAA/MTnn3+uXr16acOGDZKkxx57TNOmTVOjRo0MTgb8PDR8AAAAwE+0aNFC69ev14wZMxQUFKTly5fLYrEoIyODtg8uhYYPAAAAuImTJ0/qr3/9q1auXClJeuihhzRjxgw1b97c4GTArdHwAQAAADfRuHFjLV++XPPmzVNoaKjWr1+v6OhoTZ48WXa73eh4wE3R8AEAAAC36auvvtJLL72khQsXSrqyiXtGRoZat25tcDLg+mj4AAAAgNtUr149vfPOO1q6dKkaNmyobdu2qW3btho1apTKysqMjgdcg4YPAAAAqID8/Hy9/PLLslqtkqT4+HhZrVa1bdvW4GTAf9HwAQAAABVQp04dZWRkaM2aNWratKn27dunDh066LXXXlNJSYnR8QBJNHwAAADAHbt48aKGDh2qlJQUORwOtWnTRlarVZ07dzY6Gmo4Gj4AAADgDgUEBGjy5MnKyclRy5Yt9e9//1tdu3bVgAEDdOnSJaPjoQZj4AMAAAAqyb333qvc3Fy98sor8vDwUHJysmJjY7Vx40ajo6GG4pZOAAAAoArs2bNHiYmJysvLkyQ9//zzGjt2rIKDgw1OhpqEhg8AAACoAu3atdPu3bv1r3/9S97e3poxY4aioqK0cuVKo6OhBqHhAwAAAKrYoUOHlJiYqF27dkmS/vznPys5OVlhYWEGJ4O7o+EDAAAAqlhUVJS2bdum8ePHy8/PT/PmzZPFYtGiRYtE/4KqRMMHAAAAVKOjR4+qV69e2rRpkyTp8ccfV2pqqho2bGhsMLglGj4AAACgGkVGRio7O1vTp09XYGCgli1bJovFoszMTNo+VDoaPgAAAMAgJ06c0IsvvqhVq1ZJkh5++GGlpaWpWbNmxgaD26DhAwAAAAwSERGhlStXas6cOQoJCdHatWsVHR2tlJQU2e12o+PBDdDwAQAAAE7g3Llz+r//+z8tXrxY0pVN3DMyMtSyZUuDk8GV0fABAAAATqB+/fpatGiRlixZogYNGmjLli2KjY3VmDFjVF5ebnQ8uCgaPgAAAMDJ5Ofna+DAgZo1a5akK5u4W61WxcbGGhsMLoeGDwAAAHAyderUUWZmplavXq0mTZpoz549ateunYYPH66SkhKj48GF0PABAAAATuy7777TkCFDNHXqVElXNnHPyMhQp06dDE4GV0DDBwAAADixwMBApaSkaPPmzbr77rt16NAhde3aVX/72990+fJlo+PByTHwAQAAAC7gF7/4hfbv369BgwZJkiZMmKDY2Fht3rzZ4GRwZtzSCQAAALiY3bt3KykpSQcOHJAkvfjiixozZoyCgoIMTgZnQ8MHAAAAuJgOHTroo48+0ogRI+Tt7a3p06crKipKq1atMjoanAwNHwAAAODCDh48qMTERO3evVuS9Mwzz2jixIkKDQ01OBmcAQ0fAAAA4MKio6O1bds2jRs3Tr6+vpozZ44sFouysrKMjgYnQMMHAAAAuIlPP/1UPXv2VE5OjiTp97//vaZOnaoGDRoYnAxGoeEDAAAA3MTdd9+tjRs3KjU1VQEBAXr33XdlsVg0e/Zs0fPUTDR8AAAAgBv68ssv9cILL2j16tWSpEceeUQzZsxQkyZNDE6G6kTDBwAAALihJk2aaNWqVZo9e7bq1KmjNWvWKCoqSqmpqbLb7UbHQzWh4QMAAADc3NmzZ9W3b18tWbJE0pVN3NPT03X33XcbnAxVjYYPAAAAcHMNGjRQVlaWsrKyVL9+feXk5Cg2Nlbjxo1TeXm50fFQhWj4AAAAgBrk22+/1YABA/T2229Lktq3by+r1aqYmBiDk6Eq0PABAAAANUhISIhmz56tVatWqXHjxvroo4/Url07jRgxQqWlpUbHQyWj4QMAAABqqAsXLujVV1/VtGnTJF3ZxN1qtapDhw4GJ0NloeEDAAAAaqigoCClpqZq06ZNioyM1MGDB9W5c2cNGjRIRUVFRsdDJWDgAwAAAGq4+++/X/v379ff/vY3SdJbb72l2NhY5eTkGJwMd4pbOgEAAABctWvXLiUmJurQoUOSpL/+9a8aM2aMAgMDDU6GiqDhAwAAAHBVx44dtWfPHg0fPlxeXl6aNm2aoqKitHr1aqOjoQJo+AAAAABcV15enhITE7Vnzx5J0l/+8hdNmDBBISEhBifD7aLhAwAAAHBdsbGx2rFjh8aMGSNfX1/Nnj1bFotFS5YsMToabhMNHwAAAIBb+uSTT5SUlKQtW7ZIkrp3766UlBTVr1/f4GS4GRo+AAAAALfUsmVLbd68WSkpKfL391dWVpYsFovmzJkjOiTnRcMHAAAA4Gc5fvy4nn/+ea1du1aS9Otf/1ppaWmKiIgwOBl+ioYPAAAAwM/StGlTrV69WpmZmTKbzfrggw8UFRWl6dOny263Gx0PP0DDBwAAAKDCzpw5oz59+mjp0qWSrmzinp6ersjISIOTQaLhAwAAAHAHGjZsqCVLlmjRokWqV6+eNm/erNjYWI0fP142m83oeDUeDR8AAACASnH+/Hn1799fc+fOlXRlE3er1aqoqCiDk9VcNHwAAAAAKkVoaKjmzJmjlStXKjw8XLt27VJ8fLz+9a9/qbS01Oh4NRINHwAAAIBKV1hYqFdeeUVpaWmSpJiYGFmtVrVv397gZDULDR8AAACAShccHKzp06drw4YNatGihQ4cOKBOnTrplVdeUVFRkdHxagwaPgAAAABV6tKlS/r73/+u5ORkORwOtWzZUunp6brvvvuMjub2GPgAAAAAVIsdO3YoKSlJhw8fliT16dNHo0aNUmBgYKWsfyL/sg6cLNSBU4U6dPqCCotKVWpzyMfTpGA/H0U1ClJMeLBiGgcrok7tSrmms2PgAwAAAFBtSkpKNHLkSI0ePVrl5eVq0qSJZs6cqYcffrhC6xWX2bQi77Tm7jiu/ScLb/u8uMbBeqZzMz0a21C+3p4VurYrYOADAAAAUO3279+vxMRE7d27V5L03HPPacKECapTp85tnW+zO2TdekxTNx5VQVFZhXOY/bzV58FIJd7TXJ4epgqv46wY+AAAAAAYory8XG+99ZZGjBihkpISNWjQQKmpqXriiSduet7Rry5qUNZ+7TtRUGlZ4iPMGtc9TpH1AiptTWfAwAcAAADAUB9//LGSkpK0detWSdKTTz6pKVOmqF69etd8dtm+Uxq8JE+lNnul5/Dx9NDYP8Tq8fjwSl/bKGzLAAAAAMBQrVq1Uk5OjiZPnix/f38tWrRIFotF8+bN0w/7qfm7vlT/RblVMuxJUqnNrgGLcjV/15dVsr4RaPgAAAAAOI0vvvhCvXr10vr16yVJ3bp10/Tp0/XR1yb1X5RbLRlMkiY+2dYtmj4GPgAAAABOxeFwKDMzUwMHDlRhYaHMTVor5Klxsqn6Hqri4+mhVS/d5/Lf6WPgAwAAAOCUTp8+rb/27qNdQfeqVnjrCq/jcNhV8uVBlZz5RKVnj6rkzKeyFZ6TJIU80luB8b+57nkJTcxa/EJXl356p5fRAQAAAADgeho1aqRuA9/S/g+OyOFwyGSq2ODlKLmscwuG/uzz9n5ZIOvWY+p1X4sKXdcZ8NAWAAAAAE6puMym1E2fSVKFh73/MHn7qlbjKAV2+J3CHhskT//b2+9v6sajKi6z3dG1jUTDBwAAAMAprcg7fUebqv+HqZa/IgYuksn0374rf9Ps2zq3oKhMK/POqHu7xnecwwg0fAAAAACc0twdx3/0On/TLB0f/ahOTHpKtov51z3n3MLhOj76UZ3J7CeHrVzSlXbwh8PezzXnJzlcCQMfAAAAAKdzIv+y9p8s/NF75vuelk/9u2QvuqBvViVfc86FPStUfGyvTF61FPbbl2XyrJwbGvefLNCJ/MuVslZ1Y+ADAAAA4HQO/GTYkySTp9eVQc6rloo/36Pv9r5/9VjZ+ZMq2DhLkmR+4Dl5h0VUap6Dp67N4woY+AAAAAA4nQM3GLC8wyJkfuA5SVL+BqvKzp+Uw27TNyvGy1FeIt9m8Qps92i15XF2DHwAAAAAnM6h0xdueCyw3aPybZ4gR3mJvlkxXgUfzlXp2U/l4Ruo0G797/iJntdz8NSN8zgzBj4AAAAATqewqPSGx0wmk0K79ZeHX5BKz36qC9sXS7qyibpXYGgV5bnzp4UagYEPAAAAgNMptTluetwrIETm+5+9+rp263vl3+a+Ksxjr7K1qxIDHwAAAACn4+N589syHXabLh3ccPV16bnPZS8trsI8rjk6uWZqAAAAAG4t2M/npscv7MhSycnDMtXyl2dQXZXnn1b+hvQqzONdZWtXJQY+AAAAAE4nqlHQDY+VnD2qgi0LJEkh/+8FhXUbIMmki7mrVfTZ7irJEx1+4zzOjIEPAAAAgNOJCQ++7vv2shKdXzFesperdqt7FBD9S/k2jVVgh99Jks6vmizb5crfQuFGeZwdAx8AAAAApxPT+PoDVsGmWSo7f0KeASEK+Z8+V9+vc/9f5B3WRLZL+Tq/OuWa8+zFl2S7XHj1j3TloTCO0uIfve8ov/7TOKNddOAzORyOmz/+BgAAAAAM8LupW7T/5H/buqJj+/TVwuGSHKr35D/l16Ldjz5feu5znXl7oGQrV+hv+isg9qGrx87Oe1UlJw7e8po/PU+S4hqb9V6fe+7slzEIDR8AAAAAp/R056ZXf7YVX9T5VcmSHApI6HbNsCdJPvVbyHzvU5Kkb7NnqLzgXKXkeOYHOVwNDR8AAAAAp1RcZlPnUdkqMHDTc7Oft3YM+ZV8vT0Ny3AnaPgAAAAAOCVfb0/1eTDS0Ax9Hox02WFPYuADAAAA4MQS72mu+AizIddOaGJW4j3NDbl2ZWHgAwAAAOC0PD1MGtc9Tj6e1Tu6+Hh6aFz3OHl6mKr1upWNgQ8AAACAU4usF6Cxf4hVdY1eJklj/xCru+oGVNMVqw4DHwAAAACn93h8uN54IqbKhz6TpDeeiNHj8eFVfKXqwVM6AQAAALiMZftOafCSPJXa7JW+di0vD435fazbDHsSAx8AAAAAF3P0q4salLVf+04UVNqaCU3MGtc9zi1u4/whBj4AAAAALsdmd8i69Zimbjx6R/v0mf281efBSCXe09zlH9ByPQx8AAAAAFxWcZlNK/POaM6O49p/8vYbv7jGZj3bpam6xTR06X32boWBDwAAAIBbOJF/WQdPFerAqUIdPHVBhUVlKrXZ5ePpoWA/b0WHBykmPFjR4cGKqFNwQPZVAAAA90lEQVTb6LjVgoEPAAAAANwU2zIAAAAAgJti4AMAAAAAN8XABwAAAABuioEPAAAAANwUAx8AAAAAuCkGPgAAAABwUwx8AAAAAOCmGPgAAAAAwE0x8AEAAACAm2LgAwAAAAA3xcAHAAAAAG6KgQ8AAAAA3BQDHwAAAAC4KQY+AAAAAHBTDHwAAAAA4KYY+AAAAADATTHwAQAAAICbYuADAAAAADfFwAcAAAAAboqBDwAAAADcFAMfAAAAALgpBj4AAAAAcFMMfAAAAADgphj4AAAAAMBNMfABAAAAgJti4AMAAAAAN8XABwAAAABuioEPAAAAANzU/wfWmFmPVuYrlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 302,
       "width": 446
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%config InlineBackend.figure_format = \"retina\"\n",
    "\n",
    "G = nx.Graph()\n",
    "G.add_node(\"x1\")\n",
    "G.add_node(\"x2\")\n",
    "G.add_node(\"x3\")\n",
    "G.add_node(\"x4\")\n",
    "G.add_node(\"x5\")\n",
    "\n",
    "G.add_edges_from([(\"x1\", \"x2\"), (\"x2\", \"x3\"), (\"x3\", \"x4\"), (\"x4\", \"x5\")])\n",
    "nx.draw(G, labels={node:node for node in G.nodes})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the elements in the following dictionary represent the potential functions $\\psi_C({\\bf x_C})$ for each maximal clique $C$ which, in this case, is given by any two consecutive nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "potentials = {\n",
    "    (1,2): np.array([\n",
    "        [0.7, 0.8],\n",
    "        [0.1, 0.3]\n",
    "    ]),\n",
    "    (2, 3): np.array([\n",
    "        [0.2, 0.5],\n",
    "        [0.5, 0.3]\n",
    "    ]),\n",
    "    (3, 4): np.array([\n",
    "        [0.9, 0.9],\n",
    "        [0.5, 0.7]\n",
    "    ]),\n",
    "    (4, 5): np.array([\n",
    "        [0.7, 0.1],\n",
    "        [0.7, 0.8]\n",
    "    ])\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing $p(x_3)$\n",
    "#### The Inefficient Way\n",
    "\n",
    "$$\n",
    "    p(x_j = i) = \\sum_{x_1}\\cdots\\sum_{x_j-1}\\sum_{x_j+1}\\cdots\\sum_{x_N}p(x_1, x_2, \\ldots, x_N)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.528200000000001"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = 0\n",
    "vals = [0, 1]\n",
    "for x1, x2, x3, x4, x5 in product(*[vals] * 5):\n",
    "    Z += (potentials[1,2][x1, x2] * potentials[2, 3][x2, x3] * \n",
    "          potentials[3, 4][x3, x4] * potentials[4, 5][x4, x5])\n",
    "Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "px3 = []\n",
    "for x3 in vals:\n",
    "    px3_i = 0\n",
    "    for x1, x2, x4, x5 in product(*[vals] * 4):\n",
    "        px3_i += (potentials[1, 2][x1, x2] * potentials[2, 3][x2, x3] * \n",
    "                  potentials[3, 4][x3, x4] * potentials[4, 5][x4, x5])\n",
    "    px3_i /= Z\n",
    "    px3.append(px3_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999996"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(px3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.5813226801677082, 0.41867731983229145]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "px3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The forward/backward - passing way\n",
    "\n",
    "**Forward-passing equations**  \n",
    "$\\sum_{x_{j-1}} \\psi_{j-1, j}(x_{j-1}, x_j) =$  `potentials[j - 1, j].sum(axis=0)`\n",
    "\n",
    "**Backward-passing equations**  \n",
    "$\\sum_{x_{j+1}} \\psi_{j, j+1}(x_j, x_{j+1}) =$  `potentials[j, j + 1].sum(axis=1)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.71, 0.73])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Forward-passing message\n",
    "mu_alpha = potentials[1,2].sum(axis=0, keepdims=True).T \n",
    "mu_alpha = (potentials[2,3] * mu_alpha).sum(axis=0)\n",
    "mu_alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.07, 1.45])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Backward-passing message\n",
    "mu_beta = potentials[4, 5].sum(axis=1)\n",
    "mu_beta = (potentials[3, 4] * mu_beta).sum(axis=1)\n",
    "mu_beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999998"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "px3 = mu_alpha * mu_beta / Z\n",
    "px3.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.58132268, 0.41867732])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "px3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_pass(potential_vals, limit):\n",
    "    mu_alpha = 1.0\n",
    "    for i in range(limit - 1):\n",
    "        xi, xj = i + 1, i + 2\n",
    "        mu_alpha  = (potential_vals[xi, xj] * mu_alpha).sum(axis=0, keepdims=True).T\n",
    "    return mu_alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5282"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# An alternate way to compute Z\n",
    "forward_pass(potentials, 5).sum()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
